Normalization and Batch Effect Correction
9 May 2021
1 Objective
The objective of this notebook is to apply the downstream analysis on the filtered scATAC-seq aggregated data:
- Normalization across cells and across peaks (TF-IDF) following the method proposed by Stuart & Butler et al. 2019
- Dimensionality reduction using singular value decomposition (SVD).
- Bacht correction to reduce the dependence between the different confounders using Harmony.
- Non-dimensionality reduction (UMAP)
2 Pre-processing
2.1 Load packages
library(Signac)
library(Seurat)
library(GenomicRanges)
library(future)
library(harmony)
library(EnsDb.Hsapiens.v86)
library(stringr)
library(ggpubr)
set.seed(173)2.2 Parameters
# Paths
path_to_obj <- here::here("scATAC-seq/results/R_objects/2.tonsil_aggregated_filtered.rds")
path_to_save_obj_norm <- here::here("scATAC-seq/results/R_objects/3.tonsil_aggregated_norm.rds")
path_to_save_obj_integrated <- here::here("scATAC-seq/results/R_objects/4.tonsil_aggregated_harmony.rds")
path_tmp_dir <- here::here("scATAC-seq/2-QC/5-batch_effect_correction/1-data_asses_scATAC/tmp/")
path_to_save_dimred_uncorrect <- str_c(path_tmp_dir, "batch_uncorrected_lsi.rds", sep = "")
path_to_save_dimred_correct <- str_c(path_tmp_dir, "batch_corrected_lsi.rds", sep = "")
path_to_save_confounders_df <- str_c(path_tmp_dir, "confounders_df.rds", sep = "") 2.3 Load filtered data
seurat <- readRDS(path_to_obj)
seurat## An object of class Seurat
## 90408 features across 58049 samples within 1 assay
## Active assay: ATAC (90408 features, 0 variable features)3 Visualize UMAP without batch effect correction
# Process Seurat object
seurat <- seurat %>%
RunTFIDF() %>%
FindTopFeatures(min.cutoff = "q0") %>%
RunSVD() %>%
RunUMAP(reduction = "lsi", dims = 2:40)
DepthCor(seurat)
# Visualize UMAP
confounders <- c("library_name", "sex", "age_group", "technique", "hospital")
umaps_before_integration <- purrr::map(confounders, function(x) {
p <- DimPlot(seurat, group.by = x, pt.size = 0.1)
p
})
names(umaps_before_integration) <- confounders
print("UMAP colored by GEM:")## [1] "UMAP colored by GEM:"umaps_before_integration$library_name + NoLegend()
print("UMAP colored by sex, age group, technique and hospital:")## [1] "UMAP colored by sex, age group, technique and hospital:"umaps_before_integration[2:length(umaps_before_integration)]## $sex
##
## $age_group
##
## $technique
##
## $hospital
saveRDS(seurat, path_to_save_obj_norm)4 Run and visualize Harmony’s integration
We used Harmony integration to reduce the dependence between the main categorical batches (such as technique, sex, age group and library) of the dataset.
seurat <- RunHarmony(
object = seurat,
group.by.vars = "gem_id",
reduction = "lsi",
dims = 2:40,
assay.use = "ATAC",
project.dim = FALSE
)
seurat <- RunUMAP(seurat, dims = 2:40, reduction = "harmony")
umaps_after_integration <- purrr::map(confounders, function(x) {
p <- DimPlot(seurat, group.by = x, pt.size = 0.1)
p
})
names(umaps_after_integration) <- confounders
print("UMAP colored by GEM:")## [1] "UMAP colored by GEM:"umaps_after_integration$library_name + NoLegend()
print("UMAP colored by sex, age group, technique and hospital:")## [1] "UMAP colored by sex, age group, technique and hospital:"umaps_after_integration[2:length(umaps_before_integration)]## $sex
##
## $age_group
##
## $technique
##
## $hospital
5 Spot potential doublets
5.1 Scrublet prediction
# Scrublet
DimPlot(seurat, group.by = "scrublet_predicted_doublet")
table(seurat$scrublet_predicted_doublet)##
## FALSE TRUE
## 54837 32125.2 QC metrics
qc_vars <- c(
"nCount_ATAC",
"nFeature_ATAC"
)
qc_gg <- purrr::map(qc_vars, function(x) {
p <- FeaturePlot(seurat, features = x)
p
})
qc_gg## [[1]]
##
## [[2]]
6 Save
# If it doesn't exist create temporal directory
#dir.create(path_tmp_dir, showWarnings = FALSE)
# Save integrated Seurat object
saveRDS(seurat, path_to_save_obj_integrated)
# Save PCA matrices to compute the Local Inverse Simpson Index (LISI)
confounders_df <- seurat@meta.data[, confounders]
saveRDS(confounders_df, path_to_save_confounders_df)
saveRDS(
seurat@reductions$lsi@cell.embeddings[, 2:40],
path_to_save_dimred_uncorrect
)
saveRDS(
seurat@reductions$harmony@cell.embeddings[, 2:40],
path_to_save_dimred_correct
)7 Session Information
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS/LAPACK: /Users/pauli/opt/anaconda3/envs/Tonsil_atlas/lib/libopenblasp-r0.3.10.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggpubr_0.4.0 ggplot2_3.3.2 stringr_1.4.0 EnsDb.Hsapiens.v86_2.99.0 ensembldb_2.12.1 AnnotationFilter_1.12.0 GenomicFeatures_1.40.1 AnnotationDbi_1.50.3 Biobase_2.48.0 harmony_1.0 Rcpp_1.0.5 future_1.20.1 GenomicRanges_1.40.0 GenomeInfoDb_1.24.0 IRanges_2.22.1 S4Vectors_0.26.0 BiocGenerics_0.34.0 Seurat_3.9.9.9010 Signac_1.1.0.9000 BiocStyle_2.16.1
##
## loaded via a namespace (and not attached):
## [1] reticulate_1.18 tidyselect_1.1.0 RSQLite_2.2.1 htmlwidgets_1.5.2 grid_4.0.3 BiocParallel_1.22.0 Rtsne_0.15 munsell_0.5.0 codetools_0.2-17 ica_1.0-2 miniUI_0.1.1.1 withr_2.3.0 colorspace_2.0-0 OrganismDbi_1.30.0 knitr_1.30 rstudioapi_0.12 ROCR_1.0-11 ggsignif_0.6.0 tensor_1.5 listenv_0.8.0 labeling_0.4.2 GenomeInfoDbData_1.2.3 polyclip_1.10-0 farver_2.0.3 bit64_4.0.5 rprojroot_2.0.2 parallelly_1.21.0 vctrs_0.3.4 generics_0.1.0 xfun_0.18 biovizBase_1.36.0 BiocFileCache_1.12.1 lsa_0.73.2 ggseqlogo_0.1 R6_2.5.0 rsvd_1.0.3 bitops_1.0-6 spatstat.utils_1.17-0 reshape_0.8.8 DelayedArray_0.14.0 assertthat_0.2.1 promises_1.1.1
## [43] scales_1.1.1 nnet_7.3-14 gtable_0.3.0 globals_0.13.1 goftest_1.2-2 ggbio_1.36.0 rlang_0.4.8 RcppRoll_0.3.0 splines_4.0.3 rstatix_0.6.0 rtracklayer_1.48.0 lazyeval_0.2.2 dichromat_2.0-0 broom_0.7.2 checkmate_2.0.0 BiocManager_1.30.10 yaml_2.2.1 reshape2_1.4.4 abind_1.4-5 backports_1.2.0 httpuv_1.5.4 Hmisc_4.4-1 RBGL_1.64.0 tools_4.0.3 bookdown_0.21 ellipsis_0.3.1 RColorBrewer_1.1-2 ggridges_0.5.2 plyr_1.8.6 base64enc_0.1-3 progress_1.2.2 zlibbioc_1.34.0 purrr_0.3.4 RCurl_1.98-1.2 prettyunits_1.1.1 rpart_4.1-15 openssl_1.4.3 deldir_0.2-3 pbapply_1.4-3 cowplot_1.1.0 zoo_1.8-8 haven_2.3.1
## [85] SummarizedExperiment_1.18.1 ggrepel_0.8.2 cluster_2.1.0 here_1.0.1 magrittr_1.5 RSpectra_0.16-0 data.table_1.13.2 openxlsx_4.2.3 lmtest_0.9-38 RANN_2.6.1 SnowballC_0.7.0 ProtGenerics_1.20.0 fitdistrplus_1.1-1 matrixStats_0.57.0 hms_0.5.3 patchwork_1.1.0 mime_0.9 evaluate_0.14 xtable_1.8-4 XML_3.99-0.3 rio_0.5.16 jpeg_0.1-8.1 readxl_1.3.1 gridExtra_2.3 compiler_4.0.3 biomaRt_2.44.4 tibble_3.0.4 KernSmooth_2.23-17 crayon_1.3.4 htmltools_0.5.0 mgcv_1.8-33 later_1.1.0.1 Formula_1.2-4 tidyr_1.1.2 DBI_1.1.0 tweenr_1.0.1 dbplyr_1.4.4 MASS_7.3-53 rappdirs_0.3.1 car_3.0-10 Matrix_1.2-18 igraph_1.2.6
## [127] forcats_0.5.0 pkgconfig_2.0.3 GenomicAlignments_1.24.0 foreign_0.8-80 plotly_4.9.2.1 xml2_1.3.2 XVector_0.28.0 VariantAnnotation_1.34.0 digest_0.6.27 sctransform_0.3.1 RcppAnnoy_0.0.16 graph_1.66.0 spatstat.data_1.4-3 Biostrings_2.56.0 cellranger_1.1.0 rmarkdown_2.5 leiden_0.3.5 fastmatch_1.1-0 htmlTable_2.1.0 uwot_0.1.8.9001 curl_4.3 shiny_1.5.0 Rsamtools_2.4.0 lifecycle_0.2.0 nlme_3.1-150 jsonlite_1.7.1 carData_3.0-4 viridisLite_0.3.0 askpass_1.1 BSgenome_1.56.0 pillar_1.4.6 lattice_0.20-41 GGally_2.0.0 fastmap_1.0.1 httr_1.4.2 survival_3.2-7 glue_1.4.2 zip_2.1.1 spatstat_1.64-1 png_0.1-7 bit_4.0.4 ggforce_0.3.2
## [169] stringi_1.5.3 blob_1.2.1 latticeExtra_0.6-29 memoise_1.1.0 dplyr_1.0.2 irlba_2.3.3 future.apply_1.6.0